NOTE: This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
In the previous practice you were introduced to the concept of fields/spatially continuous data.
For this practice you will need the following:
The file contains a dataframe with geocoded observations of a series of variables. The dataframe is in tidy format, that is, each observation is a row with its corresponding variables (coordinates, variables).
The Walker Lake dataset originally was used for teaching geostatistics in Isaaks and Srivastava’s An Introduction to Geostatistics.
The dataset contains a set of false coordinates X and Y, two quantitative variables V, U, and a factor variable T. The variables are generic, but you can think of them as measurements of pollutants.
In addition, the file includes two custom functions as follows:
This is a function to obtain Voronoi polygons based on a set of points. It takes an argument sp (a SpatialPointsDataFrame) and calculates a set of Voronoi polygons. The value (output) of the function is a SpatialPolygonsDataFrame with the polygons.
This is a function to calculate k-point means. It takes a set of source coordinates (source_xy), that is, the coordinates of observations to be used for interpolation; a variable z to interpolate; a set of target coordinates (target_xy), the points to interpolate z; the number of nearest neighbors k; and a logical value to indicate whether the coordinates are latitude-longitude (the default is FALSE).
In this practice, you will learn:
O’Sullivan D and Unwin D (2010) Geographic Information Analysis, 2nd Edition, Chapters 8 and 9. John Wiley & Sons: New Jersey.
As usual, it is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:
rm(list = ls())
Note that ls() lists all objects currently on the worspace.
Load the libraries you will use in this activity:
library(tidyverse)
-- Attaching packages --------------------------------- tidyverse 1.2.1 --
v ggplot2 2.2.1 v purrr 0.2.4
v tibble 1.4.2 v dplyr 0.7.4
v tidyr 0.8.0 v stringr 1.2.0
v readr 1.1.1 v forcats 0.2.0
-- Conflicts ------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(spdep)
Loading required package: sp
Loading required package: Matrix
Attaching package: <U+393C><U+3E31>Matrix<U+393C><U+3E32>
The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:
expand
Loading required package: spData
library(plotly)
Attaching package: <U+393C><U+3E31>plotly<U+393C><U+3E32>
The following object is masked from <U+393C><U+3E31>package:ggplot2<U+393C><U+3E32>:
last_plot
The following object is masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:
filter
The following object is masked from <U+393C><U+3E31>package:graphics<U+393C><U+3E32>:
layout
library(deldir)
deldir 0.1-14
library(spatstat)
Loading required package: spatstat.data
Loading required package: nlme
Attaching package: <U+393C><U+3E31>nlme<U+393C><U+3E32>
The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:
collapse
Loading required package: rpart
spatstat 1.55-0 (nickname: <U+393C><U+3E31>Stunned Mullet<U+393C><U+3E32>)
For an introduction to spatstat, type <U+393C><U+3E31>beginner<U+393C><U+3E32>
Begin by loading the data file:
load("Walker_Lake.RData")
You can verify the contents of the dataframe:
summary(Walker_Lake.t)
ID X Y V
Length:470 Min. : 8.0 Min. : 8.0 Min. : 0.0
Class :character 1st Qu.: 51.0 1st Qu.: 80.0 1st Qu.: 182.0
Mode :character Median : 89.0 Median :139.5 Median : 425.2
Mean :111.1 Mean :141.3 Mean : 435.4
3rd Qu.:170.0 3rd Qu.:208.0 3rd Qu.: 644.4
Max. :251.0 Max. :291.0 Max. :1528.1
U T
Min. : 0.00 1: 45
1st Qu.: 83.95 2:425
Median : 335.00
Mean : 613.27
3rd Qu.: 883.20
Max. :5190.10
NA's :195
You can also check that the two functions have been loaded to your workspace.
In the preceding practice, you were introduced to three methods for interpolating spatially continuous data based on a sample of observations, namely Voronoi polygons, IDW, and k-point means.
These algorithms accomplish one of the objectives we typically have for spatially continuous data, since they provide point estimates. Recall that a field that was the outcome of a purely spatial process was defined as: \[ z_i = f(x_i, y_i) + \epsilon_i \]
A prediction of the field at a location where it was not measured was defined as a function of the estimated process and some random residual as follows: \[ \hat{z}_p = \hat{f}(x_p, y_p) + \hat{\epsilon}_p \] The first part of the prediction (\(\hat{f}(x_p, y_p)\)) is the point estimate of the prediction, whereas the second part (\(\hat{\epsilon}_p\)) is the random element.
The three methods seen previously do not, unfortunately, provide an estimate for the random element, so it is not possible to assess the uncertainty directly, which would depend on the distribution of the random term.
There are different ways in which some crude assessment of uncertainty could be attached to the point estimates. For example, a very simple approach could be to use the sample variance to calculate intervals of confidence. This could be done as follows.
We know that the sample variance describes the inherent variability in the distribution of values of a variable in a sample.
Consider for instance the distribution of the variable in the Walker Lake data:
ggplot(data = Walker_Lake.t, aes(V)) + geom_freqpoly()
Clearly, there are no values of the variable less than zero, and values in excess of 1,500 are rare.
The standard deviation of the sample is:
sd(Walker_Lake.t$V)
[1] 301.1554
The standard deviation is the average deviation from the mean. We could use this value to say that typical deviations from our point estimates are a function of this standard deviation (to what extent, it depends on the underlying distribution).
A problem with using this approach is that the distribution of the variable is not normal, and the distribution of \(\hat{\epsilon}_p\); the standard deviation is centered on the mean (meaning that it is a poor estimate for observations away from the mean); and in any case the standard deviation of the sample is too large for local point estimates if there is spatial pattern (since we know that the local mean will vary systematically).
There are other approaches to deal with non-normal variables, for instance Wilcox’s test, but some of the other limitations remain.
wilcox.test(Walker_Lake.t$V, conf.int = TRUE, conf.level = 0.95)
Wilcoxon signed rank test with continuity correction
data: Walker_Lake.t$V
V = 100580, p-value < 2.2e-16
alternative hypothesis: true location is not equal to 0
95 percent confidence interval:
418.7 475.6
sample estimates:
(pseudo)median
447.2
As an alternative, the local standard deviation could be used.
Consider the case of k-point means. The point estimate is based on the values of the k-nearest neighbors: \[ \hat{z}_p = \frac{\sum_{i=1}^n{w_{pi}z_i}}{\sum_{i=1}^n{w_{pi}}} \]
With: \[ w_{pi} = \bigg\{\begin{array}{ll} 1 & \text{if } i \text{ is one of } kth \text{ nearest neighbors of } p \text{ for a given }k \\ 0 & otherwise \\ \end{array} \]
The standard deviation could be calculated also based in the values of the k-nearest neighbors, meaning that it would be based on the local mean. Lets interpolate the field using the Walker Lake data. First create a target grid for interpolation, and extract the coordinates of observations:
target_xy = expand.grid(x = seq(0.5, 259.5, 2.2), y = seq(0.5, 299.5, 2.2))
source_xy = cbind(x = Walker_Lake.t$X, y = Walker_Lake.t$Y)
Interpolation using \(k=5\) neighbors:
kpoint.5 <- kpointmean(source_xy = source_xy, z = Walker_Lake.t$V, target_xy = target_xy, k = 5)
We can plot the interpolated field now. These are the interpolated values:
ggplot(data = kpoint.5, aes(x = x, y = y, fill = z)) +
geom_tile() +
scale_fill_distiller(palette = "OrRd", trans = "reverse") +
coord_equal()
In addition, we can plot the local standard deviation:
ggplot(data = kpoint.5, aes(x = x, y = y, fill = sd)) +
geom_tile() +
scale_fill_distiller(palette = "OrRd", trans = "reverse") +
coord_equal()
The local standard deviation indicates the typical deviation from the local mean. As expected, the standard deviation locally is usually lower than the standard deviation of the sample, and it tends to be larger for the tails, that is the locations where the values are rare.
The local standard deviation is a crude estimator of the uncertainty, again because we do not know the underlying distribution.
Other approaches based on bootstrapping could be implemented, but they are beyond the scope of the present discussion.
The issues assessing the level of uncertainty in the predictions with Voronoi polygons, IDW, and k-point means reflect the fact that these methods were not designed to deal explicitly with the random nature of predicting fields.
Other methods deal with this issue more naturally. We will revisit two estimation methods that were covered before, and see how they can be applied to spatial interpolation.
Trend surface analysis is a form of multivariate regression that uses the coordinates of the observations to fit a surface to the data.
Lets see how this would work with a simulated example. We will begin by simulating a set of observations, beginning with the coordinates in the square unit region:
n <- 180
df <- data.frame(x = runif(n = n, min = 0, max = 1),
y = runif(n = n, min = 0, max = 1))
Plot these locations:
ggplot(data = df, aes(x = x, y = y)) +
geom_point() +
coord_equal()
Lets now simulate a spatial process as follows:
df <- mutate(df, z = 0.5 + 0.3 * x + 0.7 * y + rnorm(n = n, mean = 0, sd = 0.1))
A 3D scatterplot can be useful to explore the data:
plot_ly(data = df, x = ~x, y = ~y, z = ~z, color = ~z) %>%
add_markers() %>%
layout(scene = list(
aspectmode = "manual", aspectratio = list(x=1, y=1, z=1)))
We can fit a trend surface to the data as follows. In this case, the trend is linear:
trend.l <- lm(formula = z ~ x + y, data = df)
summary(trend.l)
Call:
lm(formula = z ~ x + y, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.194786 -0.064316 0.005245 0.058750 0.230605
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.50496 0.01975 25.57 <2e-16 ***
x 0.32199 0.02496 12.90 <2e-16 ***
y 0.69977 0.02461 28.44 <2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 0.09631 on 177 degrees of freedom
Multiple R-squared: 0.8434, Adjusted R-squared: 0.8417
F-statistic: 476.7 on 2 and 177 DF, p-value: < 2.2e-16
Given a trend surface model, we can estimate the value of the variable \(z\) at locations where it was not measured. Typically this is done by interpolating on a fine grid that can be used for vizualization or further analysis, as shown next.
We will begin by creating a grid for interpolation. We will call the coordinates x.p and y.p. These we generate by creating a sequence of values in the domain of the data, for instance in the [0,1] interval:
x.p <- seq(from = 0.0, to = 1.0, by = 0.05)
y.p <- seq(from = 0.0, to = 1.0, by = 0.05)
For prediction, we want all combinations of x.p and y.p, so we expand these two vectors into a grid, by means of the function expand.grid:
df.p <- expand.grid(x = x.p, y = y.p)
Notice that while x.p and y.p are of vectors of size 21, the dataframe df.p contains {r}21 * 21 observations.
Once that we have the coordinates for interpolation, the predict function can be used in conjunction with the results of the estimation. When invoking the function, we indicate that we wish to obtain as well the standard errors of the fitted values (se.fit = TRUE), as well as the interval of the predictions at a 95% level of confidence:
preds <- predict(trend.l, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)
The interval of confidence of the predictions at the 95% level of confidence, in the form of the lower (lwr) and upper (upr) bounds:
summary(preds$fit)
fit lwr upr
Min. :0.5050 Min. :0.3109 Min. :0.699
1st Qu.:0.8325 1st Qu.:0.6404 1st Qu.:1.025
Median :1.0158 Median :0.8253 Median :1.206
Mean :1.0158 Mean :0.8241 Mean :1.208
3rd Qu.:1.1992 3rd Qu.:1.0072 3rd Qu.:1.391
Max. :1.5267 Max. :1.3331 Max. :1.720
These values indicate that the predictions of \(z_p\) are, with 95% of confidence, in the following interval: \[ CI_{z_p} = [z_{p_{lwr}}, z_{p_{upr}}]. \]
A convenient way to visualize the results of the analysis above, that is, to inspect the trend surface and the interval of confidence of the predictions, is by means of a 3D plot as follows.
First create matrices with the point estimates of the trend surface (z.p), and the lower and upper bounds (z.p_l, z.p_u):
z.p <- matrix(data = preds$fit[,1], nrow = 21, ncol = 21, byrow = TRUE)
z.p_l <- matrix(data = preds$fit[,2], nrow = 21, ncol = 21, byrow = TRUE)
z.p_u <- matrix(data = preds$fit[,3], nrow = 21, ncol = 21, byrow = TRUE)
The plot is created using the coordinates used for interpolation (x.p and y.p) and the matrices with the point estimates z.p and the upper and lower bounds. The type of plot in the package plotly is a surface:
trend.plot <- plot_ly(x = ~x.p, y = ~y.p, z = ~z.p,
type = "surface", colors = "YlOrRd") %>%
add_surface(x = ~x.p, y = ~y.p, z = ~z.p_l,
opacity = 0.5, showscale = FALSE) %>%
add_surface(x = ~x.p, y = ~y.p, z = ~z.p_u,
opacity = 0.5, showscale = FALSE) %>%
layout(scene = list(
aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
trend.plot #%>% add_markers(data = df, x = ~x, y = ~y, z = ~z, color = ~residual)
In this way, we have not only an estimate of the underlying field, but also a measure of uncertainty for our predictions, since our estimated values are bound, with 95% confidence, between the lower and upper surfaces.
It is important to note that, although the confidence interval provides a measure of uncertainty, it does not provide an estimate of the prediction error \(\hat{\epsilon}_p\). This quantity cannot be calculated directly, because we do not know the true value of the field at location \(p\). We will revisit this point later.
For the time being, lets apply trend surface analysis to the Walker Lake data.
We will first calculate the polynomial terms of the coordinates, for instance to the 3rd degree (this can be done to any arbitrary degree, however keeping in mind the caveats discussed previously with respect to trend surface analysis):
Walker_Lake.t <- mutate(Walker_Lake.t,
X3 = X^3, X2Y = X^2 * Y, X2 = X^2,
XY = X * Y,
Y2 = Y^2, XY2 = X * Y^2, Y3 = Y^3)
We can proceed to estimate the following models.
Linear trend surface model:
WL.trend1 <- lm(formula = V ~ X + Y, data = Walker_Lake.t)
summary(WL.trend1)
Call:
lm(formula = V ~ X + Y, data = Walker_Lake.t)
Residuals:
Min 1Q Median 3Q Max
-576.05 -241.79 -4.77 201.48 1055.98
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 589.5289 34.5547 17.061 < 2e-16 ***
X -1.0082 0.1903 -5.297 1.82e-07 ***
Y -0.2980 0.1727 -1.726 0.0851 .
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 292.1 on 467 degrees of freedom
Multiple R-squared: 0.06338, Adjusted R-squared: 0.05937
F-statistic: 15.8 on 2 and 467 DF, p-value: 2.29e-07
Quadratic trend surface model:
WL.trend2 <- lm(formula = V ~ X2 + X + XY + Y + Y2, data = Walker_Lake.t)
summary(WL.trend2)
Call:
lm(formula = V ~ X2 + X + XY + Y + Y2, data = Walker_Lake.t)
Residuals:
Min 1Q Median 3Q Max
-579.37 -220.93 3.96 200.66 997.30
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 320.484250 70.904832 4.520 7.86e-06 ***
X2 0.001145 0.003191 0.359 0.71994
X -0.325993 0.910872 -0.358 0.72059
XY -0.006281 0.002244 -2.799 0.00534 **
Y 3.737955 0.722097 5.177 3.37e-07 ***
Y2 -0.011409 0.002188 -5.215 2.77e-07 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 283.1 on 464 degrees of freedom
Multiple R-squared: 0.1258, Adjusted R-squared: 0.1164
F-statistic: 13.35 on 5 and 464 DF, p-value: 3.571e-12
Cubic trend surface model:
WL.trend3 <- lm(formula = V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3,
data = Walker_Lake.t)
summary(WL.trend3)
Call:
lm(formula = V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3,
data = Walker_Lake.t)
Residuals:
Min 1Q Median 3Q Max
-564.19 -197.41 7.91 194.25 929.72
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.620e+00 1.227e+02 -0.070 0.944035
X3 1.533e-04 4.806e-05 3.190 0.001522 **
X2Y 6.139e-05 3.909e-05 1.570 0.117000
X2 -6.651e-02 1.838e-02 -3.618 0.000330 ***
X 9.172e+00 2.386e+00 3.844 0.000138 ***
XY -4.420e-02 1.430e-02 -3.092 0.002110 **
Y 4.794e+00 2.040e+00 2.350 0.019220 *
Y2 -1.806e-03 1.327e-02 -0.136 0.891822
XY2 7.679e-05 2.956e-05 2.598 0.009669 **
Y3 -4.170e-05 2.819e-05 -1.479 0.139759
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 276.7 on 460 degrees of freedom
Multiple R-squared: 0.1719, Adjusted R-squared: 0.1557
F-statistic: 10.61 on 9 and 460 DF, p-value: 5.381e-15
Inspection of the results of the three models above suggests that the cubic trend surface provides the best fit, with the highest adjusted coefficient of determination, even if the value is relatively low at approximately 0.16. Also, the cubic trend yields the smallest standard error, which implies that the intervals of confidence are tighter, and hence the degree of uncertainty is smaller.
Lets compare two of these models to see how well they fit the data.
First, we use create an interpolation grid. Summarize the information to ascertain the domain of the data:
summary(Walker_Lake.t[,2:3])
X Y
Min. : 8.0 Min. : 8.0
1st Qu.: 51.0 1st Qu.: 80.0
Median : 89.0 Median :139.5
Mean :111.1 Mean :141.3
3rd Qu.:170.0 3rd Qu.:208.0
Max. :251.0 Max. :291.0
We can see that the spatial domain is in the range of [8,251] in X, and [8,291] in Y. Based on this information, we generate the following sequence that then is expanded into a grid for prediction:
X.p <- seq(from = 0.0, to = 255.0, by = 2.5)
Y.p <- seq(from = 0.0, to = 295.0, by = 2.5)
df.p <- expand.grid(X = X.p, Y = Y.p)
To this dataframe we add the polynomial terms:
df.p <- mutate(df.p, X3 = X^3, X2Y = X^2 * Y, X2 = X^2,
XY = X * Y,
Y2 = Y^2, XY2 = X * Y^2, Y3 = Y^3)
The interpolated quadratic surface is then obtained as:
WL.preds2 <- predict(WL.trend2, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)
Whereas the interpolated cubic surface is obtained as:
WL.preds3 <- predict(WL.trend3, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)
The predictions are transformed into matrices for ploting.
Quadratic trend surface and lower and upper bounds of the predictions:
z.p2 <- matrix(data = WL.preds2$fit[,1], nrow = 119, ncol = 103, byrow = TRUE)
z.p2_l <- matrix(data = WL.preds2$fit[,2], nrow = 119, ncol = 103, byrow = TRUE)
z.p2_u <- matrix(data = WL.preds2$fit[,3], nrow = 119, ncol = 103, byrow = TRUE)
Cubic trend surface and lower and upper bounds of the predictions:
z.p3 <- matrix(data = WL.preds3$fit[,1], nrow = 119, ncol = 103, byrow = TRUE)
z.p3_l <- matrix(data = WL.preds3$fit[,2], nrow = 119, ncol = 103, byrow = TRUE)
z.p3_u <- matrix(data = WL.preds3$fit[,3], nrow = 119, ncol = 103, byrow = TRUE)
This is the quadratic trend surface with its confidence interval of predictions:
WL.plot2 <- plot_ly(x = ~X.p, y = ~Y.p, z = ~z.p2,
type = "surface", colors = "YlOrRd") %>%
add_surface(x = ~X.p, y = ~Y.p, z = ~z.p2_l,
opacity = 0.5, showscale = FALSE) %>%
add_surface(x = ~X.p, y = ~Y.p, z = ~z.p2_u,
opacity = 0.5, showscale = FALSE) %>%
layout(scene = list(
aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
WL.plot2
And, this is the cubic trend surface with its confidence interval of predictions:
WL.plot3 <- plot_ly(x = ~X.p, y = ~Y.p, z = ~z.p3,
type = "surface", colors = "YlOrRd") %>%
add_surface(x = ~X.p, y = ~Y.p, z = ~z.p3_l,
opacity = 0.5, showscale = FALSE) %>%
add_surface(x = ~X.p, y = ~Y.p, z = ~z.p3_u,
opacity = 0.5, showscale = FALSE) %>%
layout(scene = list(
aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
WL.plot3
Alas, these models are not very reliable estimates of the underlying field. As can be seen from the plots, the confidence intervals are extremely wide, and in both cases include negative numbers in the lower bound. The uncertainty associated with these predictions is quite substantial.
Another question, however, is whether the point estimates are accurate. Lets add the observations to the plot:
WL.plot3 %>%
add_markers(data = Walker_Lake.t, x = ~X, y = ~Y, z = ~V,
color = ~V, opacity = 0.7, showlegend = FALSE)
Alas, the trend surface does a mediocre job with the point estimates as well.
A possible reason for this is that the model failed to capture the spatial variability. To explore this, we will plot the residuals of the model, after labeling them as “positive” or “negative”:
Walker_Lake.t$residual3 <- ifelse(WL.trend3$residuals > 0, "Positive", "Negative")
ggplot(data = Walker_Lake.t,
aes(x = X, y = Y, color = residual3)) +
geom_point() +
coord_equal()
Visual inspection of the distribution of the residuals strongly suggests that they are not random. We can check this by means of Moran’s I coefficient, if we create a list of spatial weights as follows:
WL.listw <- nb2listw(knn2nb(knearneigh(as.matrix(Walker_Lake.t[,2:3]), k = 5)))
The results of the autocorrelation analysis of the residuals are:
moran.test(x = WL.trend3$residuals, listw = WL.listw)
Moran I test under randomisation
data: WL.trend3$residuals
weights: WL.listw
Moran I statistic standard deviate = 17.199, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.4633803457 -0.0021321962 0.0007325452
Given the low p-value, we fail to reject the null hypothesis, and conclude, with a high level of confidence, that the residuals are not independent. This has important implications for spatial interpolation, as we will discuss in the following practice.
Before concluding this practice, it is worthwhile to make the following distinction between accuracy and precision of the estimates.
Accuracy refers to how close the predicted values \(\hat{z}_p\) are to the true values of the field. Precision refers to how much uncertainty is associated with such predictions. Narrow intervals of confidence imply greater precision, whereas the opposite is true when the intervals of confidence are wide.
An example of these two properties is as shown in Figure 1.
Figure 1. Accuracy and precision
Panel a) in the figure represents a set of accurate points, since they are on average close to the mark. However, they are imprecise, given their variability. This is akin to a good point estimate that has wide confidence intervals.
Panel b) is a set of inaccurate and imprecise points.
Panel c) is a set of precise but inaccurate points.
Finally, Panel d) is a set of accurate and precise points.
Accuracy and precision are important criteria when assessing the quality of a predictive model.
This concludes Practice 15.